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Abstract 

We compute the length of geodesies on a Riemannian manifold by 
regular polynomial interpolation of the global solution of the eikonal 
equation related to the line element ds 2 = gijdx t dx^ of the manifold. 
Our algorithm approximates the length functional in arbitrarily strong 
Sobolev norms. Error estimates are obtained where the geometric in- 
formation is used. It is pointed out how the algorithm can be used to 
get accurate approximation of solutions of parabolic partial differential 
equations leading obvious applications to finance and physics. 



1 Introduction 

Let (M,g) is a Riemannian manifold, i.e. a differentiable ra-dimensional 
manifold with a function g, which defines for all p £ M a positive definite 
symmetric bilinear form 

g p : T p M x T p M -> R (1) 

such that for any given vector fields X,Y € X{M) the map 

g(X,Y) : M - R, p^ g(X,Y)(p) := g p (X p ,Y p ) (2) 

is differentiable. The Riemannian metric g allows to define a metric d,M on 
M via the length of curves 

d M (x,y):= inf {L( 7 )| 7 : [0,l]^M, 7 (0)=x, 7 (l)=y}, (3) 
7 dm. 

with 

m = J ^/g J(t) m),j(t))dt. (4) 
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With this definition any connected Riemannian manifold becomes a metric 
space, and it is well known that for any compact Riemannian manifold any 
two points x,y S M can be connected by a geodesic whose length is cIm(x, y). 
If V denotes the Levi-Civita connection, then a geodesic 7 is characterized 
by the equation 

V^ 7 = 0, (5) 
which becomes (in terms of the coordinates of the values of the curve 7) 

d 2 x x „ x dx^ dx v , . 

where the well-known Christoffel symbols are 

r%, = \g Kp (d„gu P + %g w - d p9lxv ) . (7) 

This is an n-dimensional nonlinear ordinary differential equation with values 
in R n which is difficult to compute numerically in general (note the quadratic 
terms). For computing the length of a geodesic it is easier to compute the 
solution of a eikonal equation of the form 

d 2 = \Y J a *A x ) d lA 3 (8) 

(boundary conditions considered later), where (x) are functions such 

that at each x £ M. n the matrix (aij(x)) is the inverse of the positive matrix 
(gij(x)) at each point x. Here f Xi := Jjj- denotes the derivative of / with 
respect to the variable x\. In general we shall write d a f, d£f or -^sf f° r the 
multivariate derivative with multiindex a = (ai,--- ,a n )- The connection 
between the length of a geodesic which is given in local coordinates as in 
(JS|), © and the length function d 2 defined by equation (jHJ) is considered in 
section 2. This way the problem of finding the length of a geodesic is reduced 
to solving a nonlinear first-order partial differential equation in some domain 
of Euclidean space. 

The computation of d 2 is still far from trivial, however. Even if the data 
gij are analytic functions, power series expansion typically lead to power 
series solutions for d 2 with small radius of convergence. Hence, the question 
is how we can approximate the function d 2 globally. Moreover, for some 
applications such as the accurate computation of diffusions we need the ap- 
proximation of d 2 in strong norms (Sobolev norms of form H s,p for possibly 
any positive real s. For that matter recall that H°' p (R n ) = L p (W 1 ) and 
that for any s G H we may define H s ' p to be the set of all tempered dis- 
tributions <j> E § such that I- S <f) is a function in L p (M. n ), where I s is the 

s 

pseudo-differential operator with symbol <r s (£) = (l + |£| 2 ) 2 , i.e. 

J S = ?-V S 5F0, (j> £ (9) 
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If denoting the Fourier transform. The goal of the present paper can then 
be formulated as follows: find for each e > and each real s,p (p > 1) an 
approximative solution q 2 to ([8]) such that 

\\d 2 - q 2 s J s>p < e. (10) 

We shall call q 2 p an H s,p approximation to d? for reasons which will become 
apparent later. Let us motivate this ambitious task by looking at a specific 
application. There are a lot of applications for computations of the length 
of a geodesic, where applications to computations in general relativity are 
only one domain. Another important example is the leading term of the 
expansion of the fundamental solution of linear parabolic solutions (with 
variable coefficients). Varadhan showed that the fundamental solution of 
the diffusion equation 

9u — I V n- &2u 1 V b 9u fill 

(where the diffusion coefficients a^- and the first order coefficients bi in (llip 
depend on the spatial variable x only) is connected to the length d of the 
geodesic with respect to the line element ds 2 = Ylij a^dxidxj (a l i being the 
inverse of a^-) via the relation 

d 2 (x,y) = limtln.p(t,x,y). (12) 

Remark 1.1. Solving equation ([8]) we can assume that the matrix- valued 
function x — ► (aij(x)) is symmetric, i.e. aij(x) = a,ji(x) for all 1 < i,j < n. 
This is because 

d 2 (x, y) = | ^jd 2 ^ 2 . = j ^Zij 2 ( a *i + a i*) ^x^xj 

(13) 

so we can always substitute the matrix by its symmetrization ^ (a^j + aji) 
without affecting the solution d 2 . 

In [5] we have seen that for C°° coefficient functions x — ► aij{x) and 
x — ► 6j(x) and if some boundedness conditions of the derivatives are satisfied 
the fundamental solution has the pointwise valid form 

p(t, x,y) = — Ur exp ( - ^g ^ + V c fc (x, j/)t fc ) , (14) 
v 27rt V 2t fc=0 / 

where the functions x — > Ck(x,y), k > are solutions of recursively defined 
linear first order equations for each y. These equations can be solved by 
methods of characteristics or approximated by regular polynomial interpo- 
lation methods outlined in [B] . In the computation of the WKB-coefficients 
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d 2 and k > the recursive relations for c^+i involve second order deriva- 
tives of c/;, and therefore implicitly derivatives of order 2k of the squared 
metric d 2 . Hence it is of great interest to compute not only d 2 but also 
its derivatives up to a given order with high accuracy. The present work 
shows how his can be accomplished. In Section 2 we recapture some facts 
about the connection of the geodesic equation ([6]) , ([7]) and equation ([8]) , and 
prove global existence, regularity and uniqueness of the latter (family) of 
equation(s) leading us to theorem 2.3. Then in Section 3 we provide further 
analysis of the family of eikonal equations which lead us to local representa- 
tions of the solution. In Section 4 we construct first a weak approximation of 
the solution (in L p sense), and then extend this to a recursive construction 
of an .fP' p -approximation. In Section 5 we provide error estimates by using 
geometric information. Section 6 points out how the method may be applied 
for accurate approximation of diffusions, and we finish with a conclusion in 
Section 7. 

2 Global existence and regularity of the squared 
Riemannian distance d 2 

We shall only sketch the connection between geodesies and the eikonal equa- 
tion ([8]). It is almost standard, and details can be found in [5] and [I]. Our 
interest here is that the eikonal equation together with careful chosen bound- 
ary conditions has a global and unique solution. We shall have two different 
arguments for uniqueness: one is via uniqueness of an associated diffusion 
and WKB-representations (or, alternatively, Varadhan's result, cf. [8]), but 
we will have the same insight from an other point of view when we look 
at local representations of the solution in the next Section. We consider 
Riemannian manifolds where any two points can be connected by a minimal 
geodesic. For our purposes it is sufficient to consider manifolds which are 
geodesically complete. Recall that a Riemannian manifold M is geodesically 
complete if for all p € M the exponential map exp p : T p M — > M is defined 
globally on T p M. Here, T p M denotes the tangential space of the manifold 
M at p E M. The Hopf-Rinow theorem provides conditions for Riemannian 
manifolds to be geodesically complete. Especially we have 

Theorem 2.1. For a Riemannian manifold M the following statements are 
equivalent: 

• M is complete as a metric space. 

• The closed and bounded sets of M are compact. 

• M is geodesically complete. 
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Each of these equivalent statements implies that geodesies are curves of 
shortest length. Moreover, if M is geodesically complete, then any two points 
of M can be joined by a minimal geodesic. 

The connection between the arclength and equation ([5]) can be estab- 
lished as follows. First equations for minimal geodesies are obtained from 
variation of the length functional. Second Hamilton-Jacobi calculus shows 
that the length functional satisfies the eikonal equation (|8|). Since this is 
known we only sketch the main steps for convenience of the reader. Setting 
the variation of the length functional to zero we get 

L ~dh (l 2 ^) + 9ij,k ±i±j = ( 15 ) 



with L = W gij{x{r))x l xi and where we use Einstein summation. Parame- 
terizing by arclength, i.e. setting L = 1 (or r = s) we get 

2g ij x i + 2g ijt ix l x i + g ij ^x i x j = (16) 

which, upon multiplcation by g" 13 (entries of inverse of (g m j)) and rearrang- 
ing becomes the geodesic equation (HI),©. In order to show on the other 
hand that the squared length functional satisfies (JH|) we may consider the 
length functional 

l(r,x,s,y) = / L(x(u),x(u))du (17) 

Jr 

and invoke Hamilton-Jacobi calculus. This is done by introducing the vari- 
ables pi = Lj,i , and the associated Hamiltonian defined by 

H(x,p) = x % pi- L(x,x). (18) 

(here and henceforth we use Einstein summation if convenient). Then we 
may write 

x(t) = x(t; r, x, s, y) and p(t) = p(t; r, x, s, y), 
where x(r; r, x, s,y) = x and x(s; r, x, s, y) = y. and compute 

l s = -H(x(s),p(s)). (19) 

Then we may connect p to l y k by computing 



V ~ Jr \dy kPl + X dy k n * 1 dy k n P' dy k ) m 



(20) 
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(21) 



and a similar equation with respect to the variables x. Then we get the 
equations for I 2 and d 2 , i. e. the equations (jH]) and (f2"3"j) below. 

Recall that a minimal geodesic is a global distance minimizing geodesic. 
This minimal geodesic which connects x and y characterizes the Riemannian 
distance d(x, y) in an obvious way. Moreover smoothness of (x, y) — > cZ(x, y) 
for smooth diffusion and drift coefficients a« , frj follows from the following 
fact about ordinary differential equations. 

Theorem 2.2. Let F : W 1 x R n — > R &e a smooth map. Consider the 
differential system 



where x is a map JcR-* R n . TTien /or each point (xoiZ/o) ^ere exists a 
neighborhood U x V of this point and e > suc/i ifoai /or (x, w) € U X 1/ 
equation (2.69) has a unique solution x v :] — e, e[— > R n loii/i initial conditions 
x v (0) = x and x' v (0) = v. Moreover, the map X : UxVx] — e, e[— > R n defined 
by (t, x, v) — ► X(i, x, v) := £„(£) is smooth. 

Finally we get 

Theorem 2.3. Let f] C 1" some domain. The function d 2 : Q, x Q C 
R n x R n — ► R_)_ fi/ie leading order term of the WKB-expansion of a parabolic 
equation with diffusion coefficients ay) is the unique function which satisfies 
the equations 




(22) 




(23) 




(24) 



for all i,t/6 R n and with the boundary condition 



d(x, y) = iff x = y for all x, y G R n . 



(25) 



Moreover, the squareroot d is the Riemannian distance induced by 



d(x, y) := inf | y 0^(7) tVo^It : [a, b] — > R n is piecewise 
smooth with 7(a) = x and 7(6) = y j>. 
The function d 2 is a C°° -function with respect to both variables. 



(26) 
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Proof. The variation of the length functional leads to the geodesic equa- 
tion. On the other hand, Hamilton-Jacobi calculus leads us to the fact that 
the squared length functional d 2 satisfies the equation ([8]). It is clear that 
the squared length functional satisfies both equations ([23]) and ([2i]) below. 
Moreover, it is clear that the squared length functional satisfies the initial 
condition ([25]) . Uniqueness is a bit more subtle. In [8] Varadhan showed 
that 

d 2 (x,y) = lim2tlnp(t,x,y), (27) 

where p is the fundamental solution of a scalar parabolic equation with 
diffusion coefficient function x — > aij(x). Since p is unique for a strictly 
parabolic equation d 2 is uniquely determined by the equation (|27p . On the 
other hand one knows that for small t > lnp has for C°° coefficients 
a representation of type ([95]) is valid (cf. [5]). Plugging this into the 
correspondend parabolic equation leads to the eikonal equation ([8]) which 
is, hence, satisfied by d 2 . Moreover we know by V and the fact that the 
squareroot of d 2 is a metric. Hence d(x, y) = if and only if x = y, and 
the same holds for d 2 . Hence, we conclude that the global solution d 2 of 
the system of equations ([5]) . ([23]) and ([23]) is unique. Moreover, from the 
preceding theorem we can conclude that the function (x,y) — > d 2 (x,y) is 
also smooth with respect to both variables. □ 



3 Further analysis of the equation for the squared 
metric d 2 

Next we observe that the local representation of the solution of the equa- 
tions (|23p . (|24|) with the boundary condition ([25D has a local representation 
which starts with the quadratic terms. This will be used in the construction 
of a global approximation. The analysis presented here gives us two other 
insights. First, a powere series ansatz leads atmost to local and not to global 
solutions. Even if there is a local power series representation of the solution 
at each point of the domain, we do not know how a global solution can be 
constructed from this information, because we do not know the location of 
the geodesic the length of which we want to compute. If we knew, then 
computing the length would be a rather trivial task. Even the derivatives of 
the length functional would be better computed from the explicit geodesic. 
However, as we mentioned the nonlinear ordinary differential equation de- 
scribing the geodesic is harder to solve in general than the eikonal equation. 
Second, we shall see from an different point of view why the boundary con- 
dition ([25D leads to uniqueness of solutions (x,y) — > d 2 (x,y) of the system 
([23]). flM), and dS]). We have 
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Corollary 3.1. The local representation d 2 satisfying the equations (1231) . 
(|24p . together with the boundary condition (j25[) is of the form 

d 2 (x, y) = E<j a^'(y)Ax i A^ + £ |a<M ^Ax" 

(28) 

+ E| 7 |=m Jot 1 " 0) M - l ^d 2 (y + #Ax, 

JTie coefficients d a (y) are uniquely determined by a recursion obtained from 
the equations (|23j) . (|24p , /n coordinates with second order normal form, 
i.e. where d 2 is • Xi{y)Ax l Ax^ mitt Aj(y),l < i < n is the spectrum of 
{a l3 {y)), the multiindex recursion is 

+ Si S|a|>l,|7|>3,a+7=/3 ^T^O^G/H ( 29 ) 
+ Si Sa>0,|5|>3,|7|>3,a+7+5-2 l =/3 ^ilfi^siV^iv)^ . 

This confirms uniqueness. (Note that there is no loss of generality if we 
choose the normal coordinates for the second order terms). In general the 
solution is not globally analytic in the sense that d 2 is not representable by 
a globally converging power series. 

Proof. A smooth solution d 2 of the eikonal equation has the representation 
d 2 {x,y) = d(y,y) + Vd{y,y) • (x - y) 

(30) 

+ E| 7 |= 2 Jo U " ey^fd^d 2 (x + 6 Ax, y)d9. 

We abbreviate R{x,y) = £| 7 |= 2 /o^ 1 - 9) l ^ L d^d 2 (x + 8Ax,y)d6. Since 
d(y, y) = we have 

d 2 (x, y) = Vd 2 (y, y).( x -y) + R( x , y) (31) 

The 'only if'-condition of the boundary condition leads to Vd 2 (y,y) = 0. 
To see this assume that \7d 2 (y,y) ^ 0. Since R(x,y) < C||A:r|| 2 there is a 
small Ax such that Vd 2 (y,y) ■ p,Ax > C||Ax|| 2 and Vd 2 (y,y) ■ (—fi)Ax < 
— C||Ax|| 2 for some p € (0,1]. Hence there exists some p such that with 
x' := y + (pp)Ax 

d 2 (x', y) = Vd 2 (y, y) ■ (pp)Ax + R{x' , y) = 0, (32) 
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contradicting one part of the boundary condition d 2 (x,y) = iff x = y. 
Next one computes that a 13 (y)Ax t Ax 3 satisfies the equation 

d 2 (x,y) = ^^a ij (y)d*.dl., (33) 

and the uniqueness of theorem 2.3. (which we established by arguing with 
uniqueness of related diffusions and Varadhan's result, in the Atiyah-Singer 
spirit of short-range analytic expansions) identifies the coefficients a 13 (y) as 
the second order terms of local representations around y. Having obtained 
this the representation f|28H is just a multivariate version of Taylor's theorem. 
Note, however, that we do not need to invoke the uniqueness of theorem 2.3. 
but just consider a recursion obtained from a power series ansatz starting 
with second order terms. However, this would complicate the matter a 
bit so we take advantage that we know the second order terms of a local 
representation by the preceding argument. Finally we have to establish 
the recursion in (|29|) . The recursion shows directly that the higher order 



coefficients dp(y) for \/3\ > 3 are uniquely determined. Moreover, it is clear 
from (|29|) that in general the convergence radius of the full power series is 
small (if not zero). Hence in general there is no globally analytic solution 
the function d 2 : ft x ft — > K. globally analytic if for each y G W 1 the Taylor 
expansion of d 2 at y G W 1 and x G W 1 ' equals d 2 globally, i.e. 

d 2 {x,y) = ^ Migj (x-yf for all x, y G M n . (34) 

a 

Invoking the implicit function theorem equation ([5]) is equivalent to 

d^^A^^, (35) 

i 

where Aj(x),l < i < n is the spectrum of the positive (dij(x)). Since 
d 2 Xi = 2dd Xi this is equivalent to 

1 = y~]\(x)d Xi d Xi (36) 

i 

The latter equation is easier but there is no Taylor expansion around y as can 
be seen in the case of constant coefficients (and hence constant eigenvalues 
A), where the solution is 



d(x,y) 



\ i=i Ai 



9 



Remark 3.2. We use equation (I35p mainly for the theoretical purposes of 
this corollary. In general it cannot be in general used for numerical purposes 
since this would imply that we have an efficient procedure to compute the 
eigenvalue functions of a space dependent matrix. Since we are looking for 
high precision in this paper, this is not possible in general. An exception is 
the case of dimension n = 2 where we have 

, , , tT(A)(x) I ( tr(A)(x)\ 2 , , . 

Ai, 2 (x)= 2 V V 2 J - det ^)( x ) ( 38 ) 

where A(x) = (dij(x)). 

Next we plug in the power series expansion 

n 

d\x, y) = Y, A oAxf + 4(y) AxP ( 39 ) 

i=l |/3|>3 

We have 

d 2 Xt = 2\} ) (y)Ax l + £ dJ(y)ftAs"-\ (40) 
where for any multiindex (3 we define 

-(A,- | , , •••/?„) else (41) 

The term /? — 2, is defined analogously. Plugging in the power series ansatz 
and using the relation A- 1 (AM = A- 1 , this leads to 

(j2 w > 3 d 2 (y)AxP) (l-SW?) 

= (E|^|> 3 ^(!/)A^)(l-EiA) 
= (E*E H >if A^)(A^ ) 2 Ax2+ 

(42) 

+ (E, E h >i S ax-) E|/3|> 3 ^(y)AA^) 

+ (E,E a f Ax«)x 

(Ei^i^h^At^I^^^a^-'Ax^ 1 ). 

10 



This leads to 



E m > 3 d 2 (y)Ax? 

+ fa E W >i f A * a ) ( A o E|„|>s dp(y)Pi&xP) (43) 



+ (Ei E« ^a*« j x 

{Z m >3^>3 ^(y)d^y) Ax^Ax^ 



Simplifying and renaming multiindices in order to collect for multiindices of 
order (3 we get 



J2 m > 3 d}(y)AxP 

(i-EiA) ( ^* E|q|>i ( a o) -^-Ax a+2l + 



+ E* E|«|>i E| 7 |> 3 S A^(y)7iAx a+ ^ 



(44) 



+ E, E« E| 5 |>3,| 7 |>3 %5 ni d 2 (y)d 2 (y)Ax^^j . 

The latter equation leads directly to (|29p . □ 

Let us draw some consequences out of our theoretical considerations. 
There is neither an explicit solution nor leads a power series ansatz to a 
global solution in general. Neither does it help to have local solutions in 
terms of power series. Such representations are not sufficient for our pur- 
poses, since we are interested in a global solution for x — > d 2 (x,y) and do 
not know the intermediate points on the corresponding geodesic in order to 
compute the global d 2 by means of local power series representations. This 
motivates our later construction of regular polynomial interpolation of d 2 as 
seemingly unavoidable. 
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4 Regular polynomial interpolation algorithm for 
the Riemannian metric and its derivatives 



For the moment let us denote again an interpolation polynomial which ap- 
proximates the squared Riemannian distance d 2 in the L p -sense on some 
bounded domain $7 by q$ and one that approximates the squared Rieman- 
nian distance d 2 in the H s,p -sense (again on O) by q 2 p . How can we check 
that a given polynomial is an approximation in either sense? The equation 
(|8"|) gives us itself a hint how an approximation q 2 of d 2 performs. In order 
to obtain the LP error of an IP approximation q\ of d 2 we may plug in the 
approximation q 2 p into the right side of equation (|S|) and subtract the left 
side, i.e. we compute 

- ^ aij(x) Qx Qx - q 0>p = r 0jP (x), (45) 

ij 1 3 

We shall see that ro jP £ 0(h 3 ) locally (with h the mesh size of the interpo- 
lation points) implies that 

\\d(x,y) - go,p|Up(o) (46) 

converges to zero as the number of interpolation points N goes to infinity in 
such a way that the mesh size of the set of interpolation points h goes to zero. 
Note that go,p denotes the squareroot of q$ . We call an approximation q 2 p 
an i? s,p -approximation if it approximates not only d 2 in the LP sense but 
can be plugged in into all the derivatives of ([8]) of order m (i.e. multivariate 
derivatives a for \a\ < m of the eikonal equation) such that in 

9" A V ( dql p dq 2 \ &* 

[lY^^T^-J -Q^J ^v)=--r a , P (47) 

the right side staisfies ro )P € 0(/i 3+m ) locally implies that 

\\d(x,y) - q 0iP \\ H s, P{n) (48) 

converges to zero as the number of interpolation points N goes to infinity 
in such a way that the mesh size of the set of interpolation points h goes 
to zero. Accordingly, we call such q^ p (q 2 ^ p ) an LP- (H s ' p ) approximation 
of the boundary value problem 0. In the next subsection we construct a 
L p -approximation and refine the construction in the following subsection in 
order to construct iiP' p -approximations. 
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4.1 Polynomial interpolation of eikonal equation in LP sense 

We may write the eikonal equation (|5J) 

d2 ( x > V) = \ E a v<%A = \ (Vd 2 ,AVd 2 ) . (49) 

Assume that A = (dj,) is constant. The solution ([8]) with the boundary 
condition d 2 (x,y) = iff x = y is 

d 2 (x, y) = (Ax, A' 1 Ax) , (50) 

where Ax = (x — y), and A~ l =: (a^) denotes the inverse of the matrix A. 
This is easily verified by observing that 

Vd 2 = 2A~ x x. (51) 

Define 

d 2 A . 1{xj) (x,y) = Y^a lm {x ] ){x l -y l )(x m -y m ), j = 0, ■ ■ ■ ,N (52) 

ml 

we get the first recursively defined approximation algorithm for the Rieman- 
nian distance based on JV+ 1 interpolation points 
Note that the squared distance is a function 

d 2 : U x n C M n x M n -> R + , (53) 

where we define R+ := {x|x > 0}. There are several ways to approximate 
the function d 2 . In order to approximate this function we approximate first 
the function x — > d 2 (x,y), then the function x — ► c^XjX 1 ) and so on up to 
x — > d 2 (x, x ). 

We start with the approximation of x — ► d 2 (x,y). First define 
Next define 

^io(^.y) = <4-i (j/) (x,?/) + c w Uf =1 (xi - yi) 2 d 2 A - 1{xi) (x,y), (55) 
and determine a real number cio such that 

dl(xi,y) = ^J2a ij (x 1 )dl x .(xi,y)dl tX .(x 1 ,y), (56) 

i.e. the eikonal equation ([8]) with respect to x and fixed parameter y is satis- 
fied at x\. Proceeding we get a series d 2 , d?, , • • • , d 2 , • • • of approximations 
of the form 

k 

dlo(x,y) = d 2 A „ 1{y) (x,y) + J2 c joK=o n ?=i( x i ~ x l) 2 d A -i^(x, y). (57) 

j'=i 
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Having determined the real numbers cio, • • • C( k -i)a we obtain the real num- 
ber Cfco by solving 

dlo(xk,y) = \^a ij (x k )dl 0Xi (x k ,y)dl 0Xj {x k ,y). (58) 

for C&0 . Continuing this procedure for N interpolation points we get a poly- 
nomial of the form 

N 

d 2 m (x,y) = d 2 A -i {y) (x,y) + ^2c j0 U J r=0 Uf =1 (xi - xl) 2 d 2 A „ 1(xJ) (x,y). (59) 

i=i 

with N real numbers Cjo obtained recursively by plugging in d 2 with one 
degree of freedom Cjo into (f58j) . 

Analogous constructions are done to approximate x — ► d 2 (x,x^) for = 
l,--- ,iV with 

N 

d 2 Nk {x,x k ) = ^^^(x.x^+^c^n^on^^-^) 2 ^^^,)^,^), (60) 

with Cj-fc computed analogously. The construction of the functions d^ , • • • , (i 2 
suffices to approximate d 2 (we do not need to synthesize these functions into 
one function, for example by a Lagrangian polynomial interpolation). Note 
that for j = 0, • • • N the function d 2 Nk satisfies the equation 

d 2 (x,x k ) = l^2 ij a ij (x)d 2 x .(x,x k )dl.(x,x k ) 

with boundary condition (61) 

d 2 (x,x k ) = iff x = x k . 
at all interpolation points x°, ■ ■ ■ x N by construction. 

Remark 4.1. Note that in the preceding construction no restrictions on the 
choice of the interpolation points are made. This does not mean that one 
may search for an optimal choice of interpolation points and improve effi- 
ciency and convergence. We are free to choose a certain set of interpolation 
points (for example Chebyshev nodes). But these are purely computational 
aspects which will be exploited elsewhere. 

Remark 4.2. Note that we have constructed an approximation of the squared 
metric d 2 . The metric d is then approximated naturally by the squareroot 
of the approximation of the squared metric, i.e. we consider the function 

x -> d Nk (x, x k ) := yJd 2 Nk (x,x k ) (62) 

to be the approximation of the metric function x — > d(x, x k ). 
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4.2 Construction of i7 ^-approximations 

We refine the construction of the preceding section by construction of an 
approximation which solves not only , (or the set of equations (|23[) , (I24p 
with boundary conditions (|25|) ). but also all multivariate derivatives of ([H]) 
up to a given order m at the interpolation points. It turns out then that these 
polynomials are ^'^-approximations for s < m. The approximation is con- 
structed recursively again. For a multiindex (5 of order = m > 3 we de- 
note the approximations of order d 2 M ^ ^ n JV or just d 2 M ^ ^ if we do not want 
to refer to the number of interpolation points N and the dimension of the 
problem n explicitly. The choice of the mesh is free again (in principle) . We 
just assume that a set {x\, • • • ,xn} of interpolation points is given. Again 
we may construct functions x — ► d 2 M ^ (x,y), x — > d 2 M ^ (x, x 1 ),..., and 

x — > d 2 M ^ Q (x,x N ). We shall construct the first function x — > d 2 M ^^(x,y) 
for arbitrary multiindex [3. The other functions can be constructed com- 
pletely analogously. We start with the L p -approximation. 

N 

d 2 N0 (x,y) = d 2 A -!( y) (x,y) + ^2c j0 U :) r=1 Uf =1 (xi - xJ) 2 d 2 A „ 1(Xj) (x,y), (63) 

j'=i 

where the numbers Cji have been determined according to section 4.1.. Next 
we define d 2 M ^^{x, x N ) for multiindices of order \j3\ = 3. Let (3°, • • • , j3 k , • • • , j3 R 
a list of multiindices of order 3. The length R of this list is dependent of 
the dimension n of course. Start with /3° = (ffi, • • • , and let 7 be an 
multiindex with I7I = 2 such that P° — 7 = 1, for some index i. Define 
(recall that x° = y) 

d lo (x,y) = d 2 m (x,y) + -^c°p (x - yf° . (64) 

Then plug d 2 ^ 0Q (x,y) into the equation 

df-^d 2 (x,y)=df-^ , (65) 



evaluate at x = x° = y and solve for the real number c^ . Then proceed 
recursively: having defined the function x — ► d 2 p ^ k _ l ^(x,y) define 

dlo k (x,y) = dJ.^jCx.y) + c^nf-^x - x l f +1 ^(x - x k f\ (66) 

where 1 = (1, 1, • • • , 1). Then plug d 2 ^ 0k (x, y) into the equation (fB^j) . evaluate 
at y, and solve for c^ . When k = N we have got the approximation 

N , 
dlo N (x,y) = d 2 m (x,y) +£<$o n K( a - ^f +1 w (x - x k f. (67) 

fc=0 ^ 
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with N + 1 real numbers c^ for < k < N determined recursively. Note 
that the function x — > d 2 ^ 0k (x,y) satisfies the equations (jHJ) and (f69l) at all 
interpolation points x°, • • • , x N . Then we take the next multiindex (5 1 from 
the list of multiindices of order 3 (i.e. \f3 l \ = 3) where we may assume 
that (3 l - 7 1 = l fe for some multiindex 7 1 with \^\ = 2 and some index 
k. An analogous construction as in the case of f5° can be done. The only 
difference is that we start with d 2 30N (x,y) instead of d 2 NG {x,y). We get an 
approximation of the form 

N 1 

fy N (z,y) = d 2 0k (x,y) +J24^l=o(^ ~ * l f +l W^ X ~ ^ 

k=o ' ' 

where the real numbers are computed recursively by plugging the current 
approximation into the equation 

df-^d* = dr^\\^ ij a i M^), (69) 

evaluating at the current interpolation point and solving for the currently 
undetermined real number . Doing this for all the multiindices of order 
3 in the list above we get the approximation 

d M(/3 3 )( x >y) : = dpN N (x,y). (70) 

Note that by construction the function x — > d 2 M ^^(x,y) satisfies the equa- 
tion (JH|) and all its first order derivative equations 

d l x d 2 = \dU^aUx)d 2 Xl d 2 x ), l<i<n, (71) 

at all interpolation points x° = y, x 1 , ■ ■ ■ , x . This completes the stage of 
construction for multiindices of order 3. Next assume that the construction 
for the approximation 

x ^ d M(f3 m )( x >y) ( 72 ) 

of order m has been completed. Then we may list the multiindices of order 
m + 1, i.e. consider a list of multiindices d ^ 1 , ■ ■ ■ ^S Rm+1 such that \5\ = 
m + 1. The procedure is then quite similar as in the stage for multiindices 
of order 3. Therefore we give a very short description. Starting with the 
multiindex 8° there is a multiindex (3 k of order m (i.e. \(3 k \ = m) such that 
5° - [i k = U for some index i. Then we get successive approximations 

dWi^v) = d 2 M{M {x,y) + E c S° n [=o(* " x l f +X ^(x ~ xr f, (73) 

r=0 
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where the real numbers cl are succesively determined by plugging in the 
function x — ► d 2 0k (x,y) into the equation 

d^d^^feaUxKdl), (74) 

\lm J 

evaluated at the interpolation point x k (Note that d@ k = d 6 °~ li ). After 
N + 1 steps we get the approximation function x — > d 2 0N (x,y). Having 
defined x — > d 2 lN (x,y) for I = 0, ■ ■ -p — 1 the next multiindex 5 r may be 
such that there is an multiindex j3 h of order m such that 5 r — f3 h = lj for 
some index i. We may then define x — > dg Pk (x, y) 

k 

d 2 SPk (x,y) = dl I{M (x,y) + Y,c r SP U^(x - x l ) SP+1 ^(x - x r f , (75) 

r=0 

and determine the constants c r SP by plugging in the function x — > d 2 rk (x,y) 
into the equation 

^V = iW£a^)<4 m ), (76) 

V Im / 

and evaluate at x k . Finally, we get the approximation of order m + 1, namely 

4/(/3 m+l) =d 2 S H m+1N {x,y). (77) 

Note that this approximation satisfies the eikonal equation ([H]) and all its 
derivatives up to order m + 1, i.e. all equations 

dy 2 = \d« ( Y,aim(x)d 2 Xl dl m j (78) 

\lm J 

with \a\ < m + 1 at all interpolation points x 1 , • • • x N . 

Remark 4.3. Note that at some stage of the construction we may have a 
multiindex 7 such that 7 — a = lj for some a and some index Then the 
terms in the crth derivative of the eikonal equation evaluated at x k that do 
not annihilate a term of form c k H k ~Q(x — x l )' y+1 (x — x k y are quite easily 
computed. For this reason the constants of the form c k are quite easily 
computed. You can see very easily this by writing the ceth derivative of the 
eikonal equation invoking symmetry a™ = a,j. We have 

d a d 2 (x,y)=d-{lZ^ l] (x)^) 

= \ Eij (>§g) §g + \ Ey (&a l3 (x)) (79) 
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If the indicated approximation is plugged into (|79[) and evaluated at x k only 
the terms \ a ioj{ x ) (^J^ - ) fj~ (evaluated for approximations d 2 k at 
interpolation point x k ) do not annihilate terms of form c k H k ~Q(x—x l )' y+1 (x— 
x k )~<. 

5 Error estimates for the regular polynomial in- 
terpolation algorithm 

We first consider error estimates for ^-approximations, and then extend 
our estimates to .£P' p -approximations. In the whole Section we consider a 
bounded domain Q, C R n and assume that the coefficient functions Q"lj £1X6 

c°°. 

5.1 Error estimates for L p approximation 

We have 

Theorem 5.1. The approximations d 2 Nk defined in (I60p are L p - approxima- 
tions of the boundary value problems of form (|6ip . i.e. LP- approximations 
for functions of form x — > d 2 (x,x k ) for p > 1. 

Proof. Let x and y be two points connected by a geodesic curve 7 given in 
local coordinates with values in M n . Let us assume also that x and y are 
interpolation points. We have no solution for the curve 7 in general, but 
there are lets say k points z° = x, z 1 ■ ■ ■ z k = y in the image of the curve 7 
with Euclidean distance less than a certain mesh size h. Clearly, 

N 

d{x,y)=Y J d{z\z l+1 ) (80) 
i=0 

Next define an approximative distance along the geodesic of form 

n 

d g (x,y)=^2d g (z i ,z i+1 ), (81) 

i=0 

where d g is the squareroot of d 2 g (z % , z l+l ) := Ylim alm ( z m ~ z m l )(. z l ~ z i )■ 
Since y is fixed d is approximated by djyo and we estimate 

d(x, y) - d N0 (x, y) = d(x, y) - d g (x, y) + d g (x, y) - d N0 (x, y) (82) 

Our analysis showed that the local approximation of d 2 by d 2 is of order 
0(h 3 ) hence the approximation of d by d g is of order O (h^ ), hence with 
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generic constant C we have for the first summand on the right hand side of 

(El 

N 

\d(x,y)-d g (x,y)\=^2\(d(z l ,z i+1 )-d g (z\z i+1 ))\<CVh (83) 

i=0 

The modulus of the first summand on the right hand side can be estimated 
by 

\d(x,y)-<P{x,y)\<CVh (84) 

Since is a compact bounded domain, the C°° coefficient functions a 13 are 
Lipschitz Only locally Lipschitz is needed). Assuming a suitable choice of 
the points on the geodesic for the second summand we get by an elementary 
argument that 

N 

\\d g (x,y)-d m (x,y)\\ LP <Y,\\d(z\z i+1 )-d g {z\z i+1 )\\ L v<ChP-\ (85) 

i=0 

□ 

5.2 Error estimates for H s,p approximation 

Theorem 5.2. The approximations d 2 M ^ \ defined in (|77p are H s ' p - ap- 
proximations of the boundary value problems of form (I6ip for s < m, i.e. 
H s,p - approximations for functions of form x — > d 2 (x,x k ) for p > 1. 

Proof. For fixed y the function x — > d 2 (x,y) and the function x — > ^M{B m ) ( x > 
both satisfy the eikonal equation and its derivatives at any interpolation 
point by construction. That means that for all interpolation points xj, 1 < 
j < N and all derivatives 7 < m we have 

d2d 2 (xj,y) = d2d 2 M{f3m) (x j ,y). (86) 

Next recall a multivariate version of Taylor's theorem 

Theorem 5.3. If f £ C°° , then for all positive integers M we have 



/(* + v) = £i«i. (JMm ' 



*|<M 

y 1 



+M E| 7 | =M V i tia ~ 9) M ~\dnf){x + 9y)d9 



(87) 



Applying this formula, we see from our construction of x — ► d M ^ \(x,y) 
that the local order of approximation of x — ► d 2 (x, y) is 0(h 3+m ). A similar 
reasoning as in the preceding Section leads to the result. Note here that 
the same reasoning holds when y is replaced by another interpolation point 
xi. □ 
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Remark 5.4. (Sharper error estimates) A little analysis shows that the local 
order of approximation is 

h m 

d 2 {x,y) - d 2 M{M (x,y) < CP—, 
where P is the number of multiindices of order m and 



C:=2max<^ ^ sup d J d 2 (x, y), ^ sup9 7 d^ (/3) (x,y) \ (89) 

[| 7 |=M z6Q 

Similarly for \(3\ < m we get 



\-y\=M X&Q \~f\=M Xen 



l ) m-\f3\ 

d?d 2 {x,y) - d^d 2 MW Jx,y) < C {m _ m , (90) 

Since we are working on a bounded domain and d 2 is C°° there is some 
bound C, but not a priori known. However bounds for C may be obtained 
from a priori estimates by inspection of the eikonal equation. It is clear that 

d2 = lE a ^<<- (91) 
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is equivalent to 



and, hence 



d2 = iE A ^<> ( 92 ) 



\dl,\<T— , (93) 

^min 

where A m i n = minj inf xg Q A(x). Further a priori estimates for the derivatives 
may be obtained from derivatives of the eikonal equation. 



6 Analytic approximations of the fundamental so- 
lution of parabolic equations 



Consider the parabolic equation 



on some domain (0, T) x Q C (0, T) x R n , and where cc — > (aij(x)) is a 
matrix- valued C°°-function with symmetric positive matrix (a^x)) for all 
x £ Q, and x — ► b(x) is also a C°°-function. Plugging in the Ansatz 

p(t,x,y) = -^exp \- d ^1 V) +jr Ck (x,y)t k ) , (95) 
v 27rt \ 2t k=o J 
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leads to the recursive equation ([8]) for d 2 . Given d 2 the first order recursive 
equation 

- \ + \Ld 2 + \ £ (°* (*) + % j ^(*, v) = o, (96) 

together with the boundary condition 

co(y,y) = -iln v /det(ay(y)) (97) 

determines Co uniquely for each y G M n . Furthermore, having computed all 
WKB-coefficient functions q up to order k, for k + 1 > 1 the coefficient 
function c^+i can be computed via the first order equation 

(k + l)c k+1 (x, y) + | ^ + 

(98) 

= 2 Sij a «i( x ) X]z=0 Tkti~5£~ + 2" Sij a ij( x ) dxtdxj + Si ^( X ) 9^ > 

with the boundary conditions 

c fc+ i(x, y) = i? fc (y, y) if x = y, (99) 

being the right side of (f9"5|) . We will show in a subsequent paper that 
equations (i96|) ([971) . and (f98l) . (f99l) can be solved or approximated to higher 
order if ([H]) is solved or approximated to higher order. We have 

Proposition 6.1. In order to compute the WKB- approximation up to order 
k a H s,p approximation of d 2 for s > 2k is sufficient. 

Proof. In each recursion step (|98p an operator of order 2 is applied to the 
previously computed WKB-coefficients. □ 



7 Conclusion and final remarks on computational 
issues 

We have established a stable algorithm for efficient computation of the length 
of geodesies as a function of two arbitrary points on C k - Riemannian mani- 
fold with minimal geodesic as well as of partial derivatives of the length func- 
tional (of principally any order) accurately. We established error estimates 
in arbitrary Sobolev norms. We showed how the algorithm can be applied in 
order to compute fundamental solutions of irreducible linear parabolic equa- 
tions. There are many obvious applications to mathematical physics and 
finance as well as to statistics, e.g. to the maximum log-likelihood method, 
to option pricing, computing transition amplitudes etc.. Finally, we remark 
that the interpolation polynomials should not be evaluated in the way they 
are constructed. Here careful implementation of Horner schemes is needed. 
But this computational issues will be considered in a subsequent paper. 
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